home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / rwvector.lha / RWVector2.1 / src / testdfft.cc < prev    next >
C/C++ Source or Header  |  1989-08-18  |  4KB  |  140 lines
  1. /*
  2.  *    Test real FFT server
  3.  *
  4.  *    Copyright (C) 1988, 1989.
  5.  *
  6.  *    Dr. Thomas Keffer
  7.  *    Rogue Wave Associates
  8.  *    P.O. Box 85341
  9.  *    Seattle WA 98145-1341
  10.  *
  11.  *    Permission to use, copy, modify, and distribute this
  12.  *    software and its documentation for any purpose and
  13.  *    without fee is hereby granted, provided that the
  14.  *    above copyright notice appear in all copies and that
  15.  *    both that copyright notice and this permission notice
  16.  *    appear in supporting documentation.
  17.  *    
  18.  *    This software is provided "as is" without any
  19.  *    expressed or implied warranty.
  20.  *
  21.  *
  22.  *    @(#)testdfft.cc    2.1    8/18/89
  23.  */
  24.  
  25.  
  26. #include "rw/DoubleFFT.h"
  27.  
  28. overload printVec;
  29.  
  30. static void
  31. printVec(ostream& s, const DComplexVec& x);
  32.  
  33. static void
  34. printVec(ostream& s, const DoubleVec& x);
  35.  
  36. main()
  37. {
  38.   // npts must be even:
  39.   const unsigned npts = 12;
  40.   cout <<"Testing DoubleFFTServer (Double Precision FFT Server)\n\n";
  41.   
  42.   DoubleFFTServer     server;
  43.  
  44.   DoubleVec a(cos(DoubleVec(npts, 0, 2.0*M_PI/npts))); // One period
  45.   DoubleVec a2(sin(DoubleVec(npts, 0, 4.0*M_PI/npts)));    // Two periods
  46.   DoubleVec asum = a + a2;
  47.  
  48.   cout <<"**************************************\n";
  49.   cout <<"a:\n";
  50.   printVec(cout,a);
  51.   cout <<"a2:\n";
  52.   printVec(cout,a2);
  53.   cout <<"asum:\n";
  54.   printVec(cout, asum);
  55.  
  56.   cout <<"**************************************\n";
  57.   cout <<"Checking transforms of pure signals.\n";
  58.   cout <<"\nTransform of a:\n";
  59.   printVec(cout,server.fourier(a)/DComplex(npts));
  60.  
  61.   cout <<"\nTransform of a2:\n";
  62.   printVec(cout, server.fourier(a2)/DComplex(npts));
  63.  
  64.   cout <<"\nTransform of asum:\n";
  65.   // asum_fourier will be a complex conjugate even sequence
  66.   DComplexVec asum_fourier = server.fourier(asum)/DComplex(npts);
  67.   printVec(cout, asum_fourier);
  68.   
  69.   cout <<"**************************************\n";
  70.  
  71.   cout <<"Checking Parseval's theorem.\n";
  72.   cout <<"\nOriginal variance: "<<variance(asum)<<NL;
  73.  
  74.   // expand to full complex series:
  75.   double var = spectralVariance(expandConjugateEven(asum_fourier));
  76.   cout <<"Spectral variance: "<<var<<NL;
  77.  
  78.   cout <<"**************************************\n";
  79.   
  80.   DoubleVec aramp(npts,0,1);
  81.   cout <<"Checking transform of linear ramp.\n";
  82.   cout <<"\nOriginal series:\n";
  83.   printVec(cout, aramp);
  84.  
  85.   DComplexVec aramp_fourier = server.fourier(aramp)/DComplex(npts);
  86.   cout <<"Its transform:\n";
  87.   printVec(cout, aramp_fourier);
  88.   cout <<NL;
  89.   cout <<"Checking Parseval's theorem for ramp.\n";
  90.   cout <<"Original variance: "<<variance(aramp)<<NL;
  91.  
  92.   var = spectralVariance(expandConjugateEven(aramp_fourier));
  93.   cout <<"Final variance: "<<var<<NL;
  94.  
  95.   cout <<"\nBack transform:\n";
  96.   printVec(cout, server.ifourier(aramp_fourier) );
  97.  
  98.   cout <<"**************************************\n";
  99.   cout <<"Checking Nyquist.\n";
  100.   DoubleVec Nyquist(npts,1.0);
  101.   Nyquist.slice(1,npts/2,2) = -1.0;
  102.   cout <<"Original sequence:\n";
  103.   printVec(cout, Nyquist);
  104.   cout <<"Its transform:\n";
  105.   DComplexVec Nyquist_fourier = server.fourier(Nyquist)/DComplex(npts);
  106.   printVec(cout, Nyquist_fourier);
  107.   cout <<"Back transform:\n";
  108.   printVec(cout, server.ifourier(Nyquist_fourier));
  109.   
  110.   cout <<"**************************************\n";
  111.   exit(0);
  112. }
  113.  
  114. // Use special functions to print vectors, to round the output
  115. // to a fixed point number, so that differences at the machine
  116. // precision won't be seen.
  117.  
  118. static void
  119. printVec(ostream& s, const DComplexVec& x)
  120. {
  121.   for(int i = 0; i<x.length(); i++){
  122.     double rl = real(x(i));
  123.     double img = imag(x(i));
  124.     if(!(i%4) && i) s<<"\n";
  125.     s << "(" << form("%8.5f", rl) << ", ";
  126.     s << form("%8.5f", img) << ")" ;
  127.   }
  128.   s<<"\n";
  129. }
  130.  
  131. static void
  132. printVec(ostream& s, const DoubleVec& x)
  133. {
  134.   for(int i = 0; i<x.length(); i++){
  135.     if(!(i%5) && i) s<<"\n";
  136.     s << form("%8.5f",x(i)) << " ";
  137.   }
  138.   s<<"\n";
  139. }
  140.